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Based on the Green’s function (GF) equation-of-motion formalism, we develop a method to expand 
the double time Green’s function into Taylor series of the parameter A in the Hamiltonian H = 

Hq + XHi. Here Hq is the exactly solvable part and Hi is regarded as the perturbation. To restore 
the analytical structure of GF, we use the continued fraction to do resummation for the obtained 
series. The problem of zero-temperature divergence is identified and remedied by the self-consistent 
series expansion. To demonstrate the implementation of this method, we carry out the weak- as well 
as the strong-coupling expansion for the Anderson impurity model to order A^. Improved result for 
the local density of states is obtained by self-consistent second-order strong-coupling expansion and 
continued fraction resummation. 


I. INTRODUCTION 

Green’s function (GF) is widely used in the study 
of quantum many-body problems in condensed matter 
physicsii It is not only a common language to describe 
the fundamental physical concepts and processes, but 
also an important tool to do quantitative calculations for 
physical observables. Among all the methods of calcu¬ 
lating a GF, expanding it into the power series of certain 
parameter A is a basic and straightforward method for 
the Hamiltonian H = Hq -|- XHi, provided that Hq is ex¬ 
actly solvable and its GF obtainable. In cases where no 
other reliable results are available for GF, such a series 
expansion (SE) provides a reference which is accurate 
in the limit of small A. Besides quantitative informa¬ 
tion, important qualitative understanding of the system 
can also be obtained by analysing the properties of GF 
series. Well-known examples are the Fermi liquid prop¬ 
erties of weakly interacting fermions^ and the simplifi¬ 
cation of theory for lattice fermions in the large spatial 
dimension limit.— 

Various GF series expansion methods have been devel¬ 
oped so far. For a weakly interacting system, the weak- 
coupling expansion applies, where the non-interacting 
part of the Hamiltonian is chosen as Hq and Hi is the in¬ 
teraction part. Using Wick’s theorem, interacting GF is 
expanded in terms of the interacting vertex and free GF. 
Standard Feynman diagram techniques facilitate the rep¬ 
resentation and calculation of the series. Various partial 
summation methods have been developed diagrammati- 
cally, including Hartree-Fock, random phase approxima¬ 
tion, and fluctuation-exchange approximation,- 

In the other limit where the interaction is strong, the 
strong-coupling expansion can be considered, if the in¬ 
teracting part of Hamiltonian is exactly solvable. In this 
case, Hq and Hi denote the interacting part and the non¬ 
interacting part of a given Hamiltonian, respectively. Al¬ 
though the conventional Wick’s theorem no longer holds, 
various diagram techniques have been developed. Wick¬ 
like theorems were established^^— for the standard basis 
operators^ (an extension of the Hubbard operatorsiS) to 


develop diagrammatic expansions for GF. A more con¬ 
venient way is the cumulant expansion methodic intro¬ 
duced by Metzner. He expands the GF in terms of the 
hopping lines and local cumulantsi^ with unrestricted 
summations. In another elegant approach, the strong¬ 
coupling problem of original fermions is transformed 
into an effective weak-coupling problem of the dual 
fermions through a Grassmann Hubbard-Stratonovich 
transformation/^^ The GF of the dual fermions and of 
the original fermions are then obtained from standard 
Feynmann diagram technique— In both the weak- and 
the strong-coupling expansion approaches, it is necessary 
to sum partial contributions in the series to infinite or¬ 
der with a resummation method to restore the analytical 
structure of GF — If done correctly, such resummation 
can significantly extend the validity range of the theory. 

In recent years, Monte Garlo (MG) sampling methods 
have been used to carry out the series expansion of GFs, 
either in the form of determinant calculatio n^^i^° or the 
direct diagram summations—! In these methods, usually 
a large number of expansion terms can be sampled and 
summed to give the GF which is reliable in broad param¬ 
eter regimes. 

In this work, we develop a method for expanding the 
double time GF into power series of a given parameter 
A, based on its equations of motion (EOM). We call this 
method EOM series expansion. Compared to the meth¬ 
ods summarized above, this approach is distinctive due 
to following features. First, it is universal in the sense 
that the formalism does not depend on the concrete form 
of Hq. For any Hamiltonian Hq whose eigen states and 
eigen energies, or GF can be obtained exactly, one can 
always expand the full GF in terms of Hi in a recursive 
fashion. Second, there is no restriction in the operators 
A and B that define the double time GF G(A|i?)^. Sin¬ 
gle particle as well as many-particle GFs with two time 
variables can be obtained in the same general formalism. 
Third, the expansion calculation only involves double 
time GFs of Hq. The complicated multiple-time integra¬ 
tions in traditional expansion method are replaced with 
operator commutator calculations here. Fourth, in the 
present framework, the issue of partial infinite summa- 
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tion and the zero temperature divergence problem in the 
unrenormalized series expansion can be dealt with by a 
standard procedure, i.e., by using the continued fraction 
resummation and the self-consistent expansion scheme. 
Fifth, EOM of the residue of a finite order expansion is 
given, providing a possible means to estimate the error of 
the truncated series and to improve the expansion result. 

In this paper, we carry out the EOM expansion of sin¬ 
gle particle GF for the Anderson impurity model. Both 
the weak-coupling expansion to order and the strong¬ 
coupling expansion to order are obtained for single 
particle GF. In this work, we put the emphasis on the 
strong-coupling expansion. By comparing the obtained 
local density of states (LDOS) with that from the nu¬ 
merical renormalization group (NRG) method, we eval¬ 
uate the effect of the bare EOM expansion and the self- 
consistent EOM expansion, supplemented with different 
resummation methods. We show that the self-consistent 
EOM expansion together with the continued fraction re¬ 
summation gives qualitatively correct results which are 
improved in several aspects over previous ones. 

This paper is organized as follows. In Sec. II, we 
present the formal formalisms, including the the dou¬ 
ble time GF EOM series expansion, the resummation 
methods using self-energy and continued fraction, and 
the self-consistent EOM expansion. In Sec. Ill, the single 
impurity Anderson model is studied by this method. The 
weak- and strong-coupling expansions are carried out to 
second order, respectively. The strong-coupling expan¬ 
sion results obtained from different resummation meth¬ 
ods are compared with NRG results. In Sec. IV, several 
issues about the method are discussed and a summary is 
given. 


II. EQUATIONS-OF-MOTION SERIES 
EXPANSION OF DOUBLE TIME GREEN’S 
FUNCTIONS 

A. Equations of Motion Series Expansion 

We start from the EOM of retarded GFs. Let us con¬ 
sider the following retarded GF defined by two operators 
A and B at two times t and t', respectively, 

G^[A{t)m')] ^ ie{t-mA{t),B{t')]^) ( 1 ) 

Here, 9{x) is the step function. 0{t) = is the 

Heisenberg operator with respect to the Hamiltonian H. 
[A, V]j_ = XY iLYX. The plus sign is for fermion-type 
GF, and the minus sign for boson-type GF, respectively. 
{O) = Tr (e~^^/Tre~^^ is the average in thermal 
equilibrium state of H. Here h = ks = I is used. In this 
paper, the target of expansion is the GF defined in Eq. CD 
with only two time variables. We focus on the equilib¬ 
rium state where the GF depends only ont — t', although 
the method can be generalized to the non-equilibrium 


case. In the equilibrium state, the Fourier transforma¬ 
tion of G’’ [A{t)\B{t')] will be denoted as G{A\B)i^+irj, 

/ OO 

G'- [A{t)\B{t')] 

-OO 

( 2 ) 

Here rj is an infinitesimal positive number to guarantee 
the convergence of integration. 

Calculating the derivative of Eq. CD with respect to t 
or t' and transforming it onto frequency axis, one easily 
obtains the EOM for the double time GF as 

ujG{A\B)^ = {[A,B]^) + G{[A,H]\B)^ 

= {[A,B]^)-G{A\[B,H])^. (3) 

On the right hand side of Eq. new operators emerge 
from the commutator [A,H] or \B,H] and the GFs de¬ 
fined by them usually involve more particles. When the 
EOM for these new GFs are written down, even higher or¬ 
der operators and corresponding GFs will be generated. 
Usually this heirarchical EOM can not close automati¬ 
cally and approximate truncations have to be introduced 
to form a closed set of algebraic equations. In this way, 
GFs will be expressed explicitly in terms of u and some 
unknown averages. Finally, these averages will be calcu¬ 
lated self-consistently from GFs through the fluctuation- 
dissipation theorem, 

1 /■“ 1 
{BA) = — lniG{A\B)^+,^^^dcu. (4) 

Being flexible and non-perturbative, the above EOM 
formalism has been widely used in the study of quan¬ 
tum many-body systems since the early works of 
Bogoliubov,— Anderson,— Hubbard^i and others. The 
applications range from Kondo physics^^ to quantum 
magnetism— However, due to the lack of a univer¬ 
sal and systematic truncation scheme, it is difficult to 
control the precision of the resulting GFs. Especially, 
the analytical structure of GF may be violated by the 
truncation. Usually, well established truncation schemes 
are obtained empirically and only apply to specific prob¬ 
lems. Here, we will employ the EOM formalism to obtain 
a systematic series expansion for the double time GFs. 

We hrst discuss the type of Hamiltonian which is ex¬ 
actly solvable in the context of EOM. For a large class 
of Hamiltonians and operators A and B, the hierarchy of 
EOM Eq.(|3]) can form a closed set of algebraic equations. 
The GFs appearing in this set can be solved exactly 
even in the thermodynamical limit. Such Hamiltonians 
include, for examples, the non-interacting Hamiltonian 
of free electrons on a lattice Hq = and the 

Hubbard model in the atomic limit Hq — U 
The hierarchical EOM naturally close for these Hamito- 
nians. For such exactly solvable Hamiltonians, the cor¬ 
responding super-operator £, defined as CO = [Hq,0], 
has certain symmetries and hence has finite dimensional 
invariant subspaces even in the thermodynamical limit. 
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Gq{A\B)i^ can be solved exactly if operators A and B 
belong to the subspace. 

In general cases where the closure of EOM is not ob¬ 
vious, EOM of GEs can still be solved exactly if every 
eigen-state \fj.) and eigen energy of Hq can be ob¬ 
tained, Ho\^) = Efj_\^). In such cases one can construct 

the standard basis operators (SBOs)^ Act/3 = |a)(,S|. Any 

operator in the Hilbert space of Hq can be expanded by 
SBOs as 

^ ~ ^ ^ fapAcxfH, 
a/3 

B = '^gapAafj. (5) 

a/3 

Here fa /3 = (q;|A|, 5) and ga /3 = {a\B\/3). The EOM of 
GEs defined by SBOs GQ{Aai 3 \A^^)^ naturally close and 
give 

GQ{Aap\A^.h = ^ (6) 

tu -|- tja ~ t/j/S 

where {Aaa) = IZq and Zq is the partition func¬ 

tion of Hq. Gq[A\B)^ is then obtained as Go(A|H)cc = 

Xtct/3 fal3gfj.uG{Aa/3 | . 

Now suppose that the full Hamiltonian is a sum of an 
exactly solvable part Hq and a perturbation Hi. We add 
a bookmarking factor A to Hi and the Hamiltonian reads 
H = Hq + XHi. We will expand G{A\B)ui into a power 
series of A and set A as unity afterwards. Eormally, the 
Taylor series expansion of GE and averages read 

G{A\B)^ 

= Gq{A\B)^ + XGi{A\B)^ + ... + X^Gn{A\B)^ 

+Tn{A\B)^, (7) 

and 

(O) = (O)o + A(0)i -I-... -I- X^{0)n + (O)^. (8) 

Here Gi{A\B)^ and {0)i are the (i7i)®-order contribu¬ 
tions to GF and average, respectively. r„(A|H)cj ~ 
0(A”+^) and (O)^ ^ 0(A"'+^) are the residues of this 
expansion up to order n. 

Expanding the GEs and averages in Eq.dJ]) and com¬ 
paring the coefficients of A* on both sides of equations, 
one gets for i > 1 

uiGi{A\B)f^ 

— ([^) B\^)i + Gi_i([A, + Gi([A, Hq\\B)i.j 

= {[A,H]±), - G,-i{A\[B,Hi])^ - G,{A\[B,Hq])^, 

(9) 

and for i = 0 

a;Go(A|H)cc = {[A, B]^)q + Go{[A, Hq]\B)i^ 

= {[A,B]^)q-Gq{A\[B,Hq])^. (10) 


The residue r„(A|i?)cj of the n-th order expansion satis¬ 
fies the EOM 

= {[A, + G„([A, Hi]\B)^ + r„([A, H]\B)^. 

= {[A,H]^)«-G„(A|[H,i7i])cc-r„(A|[H,i7])c.. 

( 11 ) 

In this formalism, one can choose to use the left-side or 
the right-side EOM formula differently for different order 
i. 

To solve the averages involved in the above equations, 
expanding Eq.([4|) gives 

{BA)i = -/ lmGi{A\B)^+ir,-^——duj (12) 

for 0 < i < n. The residue (O)^ is obtained from 

1 1 
{BA)^ = / Imr„(A|H)cc+i7) duj. (13) 

Since Hq is exactly solvable in the sense discussed 
above, the EOM of Gi{A\B)i^ in Eq.(l9l) will close because 
the series A, [A,i7o]) [[-^t Hq], Hq], ... generates closed set 
of operators. In cases where this is not obvious or too 
complicated, one could decompose A and B into SBOs 
and study the Taylor series expansion for the GEs defined 
with SBOs. In any case, Gi(A|H)(^ can be expressed in 
terms of the lower order GF Gi-i{A']B)i^ with more com¬ 
plicated operators A'. Repeatedly employing Eq.®, one 
can then reduce Gi_i(A'|i3)^ to Gi_ 2 (A"|H),^, and so on. 
Finally the GF component of every order z > 1 can be 
reduced to the type Go(A|H) with different operators A. 
These zeroth order GEs are exactly solvable. Therefore, 
Eq.([n]) provides a practical way of calculating arbitrary 
order contributions to G(A|H),^. 

The n-th order residue of the expansion rji(A|H),^ is 
determined by its EOM Eq. (ED which cannot be solved 
exactly. By truncating the hierarchical EOM Eq. ED, 
one could obtain an approximate result for r„(A|i?)^. 
This result could be used to evaluate or accelerate the 
convergence of the expansion. 

Equation (9) contains {[A,B]±)i, the z-th order con¬ 
tribution to {[A,B]±). It must be calculated through 
Eq.(ED> which leads to a set of algebraic equations for 
the averages. In the conventional GF EOM approach 
with truncation approximations, neither consistency nor 
sufficiency is guaranteed for this set of equations. In the 
rigorous series expansion presented here, for each order 
z, the set of equations is both consistent and sufficient. 
That is, a unique solution of the averages at z-th order is 
always obtainable. The n-th order residue of the average 
{[A,B]±)^ needs to be calculated self-consistently using 

Eq.(ED- 
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B. Resummation Methods 


Truncating the Taylor series of a GF at a finite order 
always produces the problem of causality. This is most 
easily demonstrated by Taylor expanding the Lehmann 
representation of G{A\B)i^ (taking the fermion-type GF 
as an example) 


G{A\B)^ 



m,n 


_|_ g-/3£„ 

W -|- Em — En 


X 


mn • 


(14) 


Here \m) and Em are the eigen state and the eigen en¬ 
ergy of E[, respectively. Z = J2m is the partition 

function, and Xmn = {rn\A\n){n\B\rn) is the matrix ele¬ 
ment. It is seen that GF has only real simple poles and 
(3 appears on the exponent. 

We can formally expand Z, Em and Xmn into power 
series of A to obtain 


z = Y,xz, 

oo 

Em=Y.yE^^ 

oo 

Xmn = Y.\^X^±. (15) 

i=0 

The Lehmann representation for the z-th order term 
Gi{A\B)^ can be obtained by putting these expansions 
into Eg. (ITT)) and collecting terms proportional to Ah The 
first two orders read 


GM\B)u. 


and 


1 


E 




^ ^ y(o) 

, , , p(0) p(0) 

m n ^ \ -C/m -C/n 


E 


(16) 


Gi{A\B\ 

1 e 




_y(0) _ ^ p(l) y(0) . y(l) 

^^x\'Ynn -^mn \ ^mn 


E 


-^E 




. , _|_ Z7'(0) 

iU + I^rn ~ -^n 

_Zi y (0) _ o mCl) y (0) , y (1) 

^0 HJ-'m -^mn “r x\mn 


t 1 -L 

U) -\- h/ra ~ -t^n 




E^B _ 


u + E^°^ - E^^'>' 




(17) 


Eg. dTBl) shows that Go(H|i?),^ has real simple poles 
which gives the expected causality of the retarded GF 
Go(H|il)^+i,,. However, Gi(H|i?),^ has real second-order 
poles and thus violates the causality. Moreover, some 
terms in Gi(H|i3),^, including Zx = -/3 


have a factor /3 which comes from the Taylor expansion of 
Em in the Boltzmann factor. Thus Gi(H|il)^ will diverge 
in the zero temperature limit unless (3 factors in different 
terms cancel. We call this problem the zero-temperature 
divergence problem, which occurs when the ground state 
of E[ is unanalytical at A = 0. In general, the fc-th order 
term Gk{A\B)^ has (fc -|- l)-th order poles and contains 
terms with factor j3^. As a result, a truncated series of 
GF almost always violates the analytical structure and 
diverges at zero temperature. In order to avoid these 
problems, it is necessary to sum part of the expansion 
contributions to infinite order. 

In the following, we discuss possible resummation 
methods to recover the analytical structure of GF. The 
zero-temperature divergence problem will be considered 
in the next subsection, invoking the self-consistent EOM 
expansion. For interacting electron systems, a conven¬ 
tional practice is to directly expand the self-energy (SE) 
up to some hnite order and then insert it into the Dyson 
equation to produce GF. Suppose we have obtained the 
series expansion of the single-particle GF G{cka\cl.^)ui for 
a lattice Hamitonian up to order A" 

G{cka\cl.^)ui 

ss Go(CfcCT |cj)^)i^ + AGi (cfco-|c|.^)aj + ...+ A GniCkalCf.^)^- 

(18) 


Here Cka and are electron creation and annihilation 
operators, respectively. The GF obtained from the SE 
resummation is denoted as GsE{,Cka\cl^)ui, 


GsE{cka\clju, 


1 

Gq (fc, w) — E(c/cct|c^o-)‘-J 


(19) 


Here, t/^^(fc,w) = oj + p, — €k is the exact GF of the 
non-interacting Hamiltonian H{V = 0). Here V rep¬ 
resents the interaction strength in H. In Eo. dT^ . we 
substitute the SE with its formal series truncated at n- 
th order E(cfe,^|c[^)^^ « ^o{cka\cl^)ij + XEi{cka\c\^)ui + 
...-I-A"S„(cfc(T|c[^)(^. GQ^(k,u:) should also be expanded 
if H{V = 0) depends on A. Expanding GsE{cka\c'^.^)u) 
into power series of A and comparing it with the series of 
G(ck(T\cl^)u! to order A", we can fix the expansion terms 
of SE. GsE{Ckcr\cl^)uj obtained in this way is exact up to 
order A" but contains approximate terms from A"“'"^ to 
A°°. Note that other equations derived from the diagra- 
matic resummation, such as the Larkin equation^ could 
also be used to do resummation in a similar way. 

In the strong-coupling expansion of the Anderson im¬ 
purity model, the SE resummation method has been used 
to calculate the impurity GFi^ The resulting GF vio¬ 
lates the causality and does not obey the sum rule (see 
below). Similar problems such as negative spectral func¬ 
tion also appear in the weak-coupling expansions ^2^ in¬ 
specting the analytical structure of SE E,{cka\c\^)ui = 
Cofc + Sm ~ <-^mk)r ^e find that a truncated se¬ 

ries of SE violates the correct analytical structure and has 
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zero-temperature divergence problem, if only the poles 
{uJmk} and weights {c^fc} contain A and /3. In special 
cases where the poles have A*-" contributions only, the 
SE expansion to order A"”^ does not produce non-simple 
poles and the causality will be fulfilled. As will be shown 
below, this is the case of the weak-coupling expansion 
to order U'^. We conclude that in general, the resum¬ 
mation method from the truncated bare expansion of SE 
can guarantee neither the correct analytical structure nor 
the correct zero temperature limit. 

To overcome the problem with analytical structure, 
Pairault et al. suggested a resummation method based 
on the continued fraction (CF).— For the single parti¬ 
cle GF G{ckcT\cl^)uj, one can construct a CF of the form 


GcF{da\dl)u 


ao(A) 


w -I- 61(A) - 


ai(A) 


w -I- 62(A) — 


a2(A) 


OJ + 63(A) - ... 


( 20 ) 


It was proven that for real parameters a; (A) > 0 and 
real 6;(A) {I = 1, 2,...), the above expression always gives 
the correct analytical structure of GF, i.e., it consists of 
real simple poles. This can be understood by regard¬ 
ing Ea. (|20l) as the local GF of a free semi-infinite chain 
Hamiltonian. To carry out the GF resummation, we for¬ 
mally expand the coefficients a;(A) and 6;(A) (Z = 1,2,...) 
into Taylor series of A up to A” order, and then compare 
the obtained GcF{da\d^rT)‘-i G{da\dl)ui to order A" 
to fix the expansion coefficients of a;(A) and 6;(A). Usu¬ 
ally, due to the finite number of poles produced by the 
GF expansion, a CF with finite levels is sufficient for this 
task. Besides being exact up to order A", the obtained 
GF Gcf (da|4)- is causal. This resummation method 
was used in the strong-coupling expansion study of two 
dimensional Hubbard modelddrii 


C. Self-Consistent EOM Series Expansion 


effect was attributed to the fact that the intra-Hubbard- 
band hopping process of electrons involves the energy 
scale t/T instead of t/U— and cannot be described accu¬ 
rately hy t/U expansion. Generally speaking, this prob¬ 
lem reflects that the ground state energy or the density 
matrix is not expanded on the same footing as the exci¬ 
tation energies. Inspecting the Lehmann representation 
of SE shows that this problem also exists in the SE re¬ 
summation. 

In the bare expansion formalism Eq. Q , the only place 
where /3* emerges is the component {[A,B]±)i {i > 1). 
Gi{A\B)i^ generally has {i -I- l)-th order poles. When 
{[A,B]±)i {i > 1) is calculated from the same order 
GF component using the fluctuation-dissipation theorem, 
/3* will be produced through the frequency integration. 
Therefore, this problem could be solved if the averages 
are not calculated order by order from each Gi{A\B)i^, 
but from the full G(A|i?)^ after its correct analytical 
structure has been restored. Below, we propose the self- 
consistent GF EOM expansion scheme in this spirit. 

We start from the EOM for the full GF Eq.([3]). For 
simplicity, here we consider the left-side EOM only. The 
corresponding formula for the right-side EOM can be de¬ 
rived similarly. We define the renormalized zeroth-order 
GF Go{A\B)^ by the EOM 

ojGo{A\B)^ = {[A,B]^) + Go{[A,Ho]\B)^. (21) 

Unlike in the bare expansion Ea. ifTOl) . ([A, H]^) here is 
the thermodynamical average with respect to full Hamil¬ 
tonian H. It will be calculated self-consistently from 
the CF-resummed GF which has correct analytical struc¬ 
ture. Subtracting Eq. m from the the left-side EOM 
Eq.Q and defining the renormalized n-th order residue 
as r„(A|H)^ = G{A\B)^ - Go{A\Bh - GiiA\B)^ - ... - 
Gn{A\B)^, we get the EOM for the zeroth-order residue 
ro(A|H)^ as 

ujro{A\B)^ = Goi[A,Hi]\B)^+Toi[A,H]\B)^. (22) 

For the next order renormalized GF Gi{A\B)i^, we re¬ 
quire that it satisfy the EOM of ro(A|i?),^ at the leading 
order of A, i.e., Ea. d^ with H replaced by Hq. We have 


The CF resummation method can overcome the causal¬ 
ity problem in the truncated series of GF and recover 
the correct analytical structure. However, it still has 
the problem of zero-temperature divergence. Using 
Lehmannn representation, we have shown that in gen¬ 
eral, P factors will appear in Gi{A\B)i^ {i > 1). The CF 
resummation procedure transmits these /3 factors into 
the parameters ai{\) and 6;(A), leading to an unphysi¬ 
cal shifting of certain poles to infinity as T approaches 
zero. Indeed, in the strong-coupling expansion study 
of a two dimensional Hubbard model^^— /3 factors ex¬ 
plicitly appear both in the bare expansion terms of GF 
and in the parameters of CF. The validity range of the 
CF-resummed GF is thus limited to temperatures higher 
than some energy scale. For the Hubbard model, this 


uJG^{A\B)^=Go[[A,H^]\B)^ + G^{[A,Ho\\B)^. (23) 

Similarly, EOM of the residue ri(A|H)^ = ro(A|i3)(^ — 
Gi{A\B)^ can be obtained by subtracting Ea. d^ from 
Eq. (1^ . This procedure is carried out repeatedly to pro¬ 
duce EOM of the renormalized GF component for f > 1, 

ujGPA\B)^ = G,-i{[A,HP\B)^ + GP[A,Hq]\B)^. (24) 

For z = 0, Eg. d^ applies. The n-th order residue 
r„(A|H)(^ satisfies the EOM 

ccr„(A|H)<, = G„([A,Hi]|H)<,+r„([A,iL]|H)„. (25) 

Series expansion using the right-hand side EOM can be 
derived similarly. Note that in this self-consistent series 
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expansion, one can use either the left-hand side or the 
right-hand side EOM formula throughout the derivation, 
but cannot mix them in different orders. This differs 
from the bare EOM expansion. 

This EOM expansion scheme only involves full av¬ 
erages ([A, which need to be calculated self- 

consistently from the corresponding full GEs. To be free 
from the zero-temperature divergence problem, these GE 
must have real simple poles and therefore CF-resummed 
GF should be used. 

A formal solution of the above self-consistent EOM ex¬ 
pansion can be obtained in terms of the Liouville operator 
C as. 




= G,_i ( A- ^AB] 

LO — Lq 


r„(A|B)^ =G„ A--A 

LO — L 


B 


(26) 


In the weak-coupling expansion, we decompose 
^Aim — Hq T 

k<7 k<7 <y 

(28) 

and 

Hi = Un^ni-'^aaTia- (29) 

(7 

Here, ida = Cd + ao-. Oo- is a parameter to be determined 
by the principle of convenience. For an example, its value 
could be fixed by requiring that the first-order contribu¬ 
tion Gi{da\dl)c^ is zero. In that case, expanding GF into 
Taylors series of Hi amounts to perturbation around the 
Hartree-Fock Hamiltonian. 

To facilitate comparison, we will apply the bare EOM 
series expansion and use the SE resummation. Here the 
GF is of fermion type and the anti-commutator of oper¬ 
ators A and B is denoted as {A,B}. The zeroth-order 
GF GoidMh is easily solved from its EOM, 


Here, Co and Ci are the Liouville operators of Hq and 
Hi, respectively. They act on any operator O as CqO = 
[6, Ho] and £i 6 = [6, Hi]. 


III. WEAK-COUPLING EXPANSION FOR 
ANDERSON IMPURITY MODEL 


In this section, we apply the formalism developed 
above to the single impurity Anderson model. This 
model is one of the best studied models for correlated 
electron systems, due to its importance in the dilute mag¬ 
netic impurity problem, in the quantum dot physics, as 
well as in the application of dynamical mean-field the¬ 
ory for Hubbard model. The Hamiltonian of the single 
impurity Anderson model reads 


Ha 


im 


fJ-)c].^Ck(T + 'y ] Vfco- A d^^Cka-'^ 

k<7 k<7 


(Cd - /r) ^ rio- + Gn-i-n4,. 

CT 


(27) 


Here we consider spin-dependent energies and hybridiza¬ 
tions of bath electrons. Ua = dJ-A is the electron 
number operator of impurity. The influence of bath 
to impurity is described by the hybridization function 
Aa{uj) — J2k ~ ^krr)- Throughout this work, we 

set the chemical potential p, = 0 as the zero point of 
frequency. 

Both the weak- and the strong-coupling expansion for 
this model have been obtained before. Here, for demon¬ 
stration purpose only, we apply the EOM series expan¬ 
sion method to derive the weak-coupling expansion for 
the local GF to order, to recover the well-known re¬ 
sults of Yamadcfc^ at the same level. 


Go(A|4)- 


1 

UJ-eda- ro-(w)' 


Here r,^(a;) = J2k ^fca/4 “ <^ka)- 

For i> 1, the EOM for Gi(A|4)i.j reads 


(30) 


ujGMdi)u. 

= + g,_i([a, Hi]i4)^ + a([a,Ho]|4)<.. 

(31) 

Using (1)*>1 = 0, [A,Ho] = J2kAkCka + idada and 
[A, Hi] = Uriada - Oct A, we get 

(w - eda)Gi{dcr\dl)u; 

= ^I4A(cfe,|4)^-a,A_i(A|4)- 

k 

-|-HGi_i (rig-A |4)‘.>j' (32) 

We use the left-side EOM for the i-th order new GF 
Gi(cfcg|4)w to obtain 


G,{cka\dih = —^A(A|4)-- (33) 

W - Ckcr 

Putting Ea. (l33l) into Eg. (1321) . we get 

a(ai4)- 

= Go(A| 4)- [HA_i(ngA|4)<. -a<.G,_i(A|4)<.] • 

(34) 

For 1 = 1, this equation involves Go(ngA|4)t.jj which is 
easily solved from its right-side EOM as Go(ngAl4)<.>j = 
(ng)oGo(Al4)<.-j- Eg.dMl) then produces 


Gi(ai4)- 


U {klfj^O (X(y 

[w - ida - rg(w)]^ 


(35) 
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For i = 2, Ea. (l34)) involves a new GF Gi{nada\d'l)u,- 
After some calculation, its EOM is solved to produce 

G\{Ti^d(y\ d^ ) q; 

= [{n^)i + UK^iu;)]Goid^\di)^ 

d~{Tl(f)o [U {Tla}o CKcr] Gq (d^r jd],. )^ . 

(36) 


Now we can choose to simplify the expression. By 
assigning = U{ng)Q^ the expansion gets simplified as 

GQ{da\d\)uj = [w - eda - rcr(a;)]“^; 

Gi{da\d]j)ui = 0; 

G2{dM)u. = Gl{dM)u. [U{nd)i + U^K^iu:)] . 

(43) 


Details of this calculation is presented in Appendix A. 
Putting Ea. (IM)l and Ea. (l36)l into Ea. (IM|) we obtain 


G2{dM)^ = U[{n^)^ + UK,{uj)]Gl{dM\ 
+ [U{nd)o - Gl{dM)u,- 

The function Ka(u}) appearing here reads 

A'cr(w) 


(37) 


dOff (ei)doff (e2)poo-(e3) 
w + Cl - £2 - £3 


■F'(£i,£2,£3)d£ide2d£3, 

(38) 


where poai^) = — l/7J'IwGo(dCT|dJ.)e+i^ is the free local 
density of states of spin a and 


The corresponding SE reads 
S,(a;) = U [{n^)o + (n<,)i] + U^K^{uj) + 0(C/3). (44) 

The SE-resummed GF is finally obtained as 

GsE{da\dl)uj = [g^J-iuj )^ Notc that 
(ng-)i has a factor jd and may lead to zero-temperature 
divergence in SE. Here, the causality problem does not 
appear since the obtained SE has real simple poles only. 
The reason, as analysed in general in Sec.IIB, is that the 
poles of SE for the Anderson impurity model contain 
terms of the order C/*-^ and expansion to order 
does not produce non-simple poles. In the paramagnetic 
and particle-hole symmetric case, {ng)^ = 1/2 and 
{ng)i = 0, the zero-temperature divergence problem 
does not appear. The above expression is equivalent to 
the Matsubara SE obtained by Yamada,— 


^’(£ 1 , £2, £ 3 ) = ^^(£3) [^(£2) - n(ei)] -f n(ei) [1 - 71(62)] , 

(39) 

with n(e) = l/(e^'^ -1- 1) being the Fermi-Dirac distri¬ 
bution function. The unknown averages involved in the 
above expressions are ( 77^)0 and ( 71 ^) 1 . From 

1 1 

I ^a)a;-|-z?7 

^ J —00 ^ 

they are calculated as 

/ oo 

/0off(e)7i(e)de; 

-00 

(no-)! — {JJ (71ct)o £^a) ^ 

f°° Po<t(£i)po<t(£ 2) r 7 ^ / M j j 

/ -[n(ei) - 71(62)] d£id£2. 

j -00 £1 ~ £2 

(41) 

Until now, we have obtained the lowest three orders 
of GF, Gi{da\d\)i^ {i = 0,1,2). To carry out the SE 
resummation, we insert the truncated expansion of SE 
Ecr(a;) Eocr(w) + AEiCT(a;) -|- A^S 2 ct(w) into the Dyson 
equation GsE{dg\dl)^ = [GoJiuj) - ^^,( 0 ;)] , expand it 

into a Taylor series of A, and compare the A* term with 
Gi{dg\di)uj (i = 0,1,2). Goai^^) = 1/[w - 6d - r,^(w)] is 
the non-interacting local GF. The SE is given as 

E0{7(^) ^(Ti 

Ei,(a;) = G^2(^^|^t)^Gi(d,|4)^; 

£ 2 .( 0 ;) = Go"(d„|4),,G2(d.l4)^ 

-G^\dg\di)^Gl{dg\di)^. (42) 


IdgiiuJr,) = ^ + U^ f G3(T)e“""dT; 

^ Jo 

GoiiiOn) = [iiOn - ro.(7a;„)]“^. (45) 

When combined with the dynamical mean-field the¬ 
ory, the above SE produces the iterative perturbation 
theory^ (IPT) which describes the Mott metal-insulator 
transition very well. For the particle-hole symmetric case 
and paramagnetic bath, the above bare expansion with 
SE resummation has neither causality problem nor zero- 
temperature divergence problem due to {n-g) = 1/2. This 
simplest form of IPT fails, however, away from particle- 
hole symmetry or in the magnetic bath. Extension of the 
original IPT to such situations received some research 
effort^ in the spirit of interpolation between various ex¬ 
act limits. Ea. (l4^ reminds us that the zero-temperature 
divergence problem may occur in the bare expansion. 
The recovery of atomic limit is also difficult by SE- 
resummation away from particle-hole symmetry. These 
problems could be remedied by the self-consistent EOM 
expansion supplemented with GF resummation. For the 
moment, we leave in depth discussions of this issue to the 
future and focus on the strong-coupling expansion of the 
GF. 


IV. STRONG-COUPLING EXPANSION FOR 
ANDERSON IMPURITY MODEL 

In this section, we carry out the strong-coupling ex¬ 
pansion for G(dcr|4)a; to the order V)?. Besides test¬ 
ing the EOM expansion methods, the obtained formula 





may also serve as a useful strong-coupling impurity solver 
for the dynamical mean-field theory to describe the anti¬ 
ferromagnetic insulating phase, for which existing weak- 
coupling-based theories such as IPT and the functional 
renormalization group method^ are faced with difficul¬ 
ties. 

We investigate the impurity spectral function, electron 
occupation, and the double occupancy at the particle- 
hole symmetric point both for the paramagnetic and the 
magnetic cases. In the paramagnetic case, our results 
agree with the direct expansion results We compare 
the local density of states from three different calculation 
schemes: bare EOM expansion supplemented with SE re¬ 
summation, bare EOM expansion with CF resummation, 
and the self-consistent expansion with CF resummation. 

To do the strong-coupling expansion, we decompose 

HAim — 1^0 “b ^1 into 


given by G(d„|4)^ = Ea/3 

The EOM of the zeroth-order GF reads 


ujGo{A'^p\A^g)i^ — {{A'^p, A'^g})o+Go{[A’^p, Hq] 

(49) 

Using the relations = SpsAaj + S^aAsp, 

(^a/3)o = 5ape~^^°‘lZQ and A'^p,HQ 

{Ep — Ea) A'^p, we obtain 

fla + 0 ,p 


Go{AIp\A:^ 1 )^ = 5^^5ps 


w -f Ea — Ep 


(50) 


Here aa = e“^^“/Zo, and is the parti¬ 

tion function of EIq. 

For f > 1, the EOM reads 


ujGi{Al^p\A''^\)^ — -I-Gi_i([A^^,Hi] 1^4), 


i?o = E' 


(^kac\^Ckcr + u n^nj, + ed ^ i 


ka 


and 


Hi — 'y ( Vka + dl^Cka^ ■ 

k (7 


(46) 

+Gi([A^^, Ho 



To simplify the notation, we 

expand 


commutators involved in \A'^p,Hi 

as 

(47) 

L 



{A^^p, 4 ,} — y^ 


(51) 


Obviously, E[q is exactly solvable. Due to the existence of 
Un-fTii in Hq, the hierarchy EOM for the local GF closes 
at the second level and it is more convenient to work with 
SBO formalism. We denote the eigen-state and eigen- 
energy of ho = Un^ni + CdEo-^o- ^.s |a) and Ea, re¬ 
spectively. That is, ho\a) = Ea\a). The SBOs {Aap} 
are defined as the projector operators Aap = |q!)(/ 3|. 
They satisfy the algebraic relations AapA^^, = Sp/^Aai, 
and Ea = 1- The latter plays the role of the kine¬ 
matic sum rule^i and is important for the self-consistent 
solution for the averages. Eq is now written as 

Hq = y ] ^ka-C/ca^ktT + ZeM(A (48) 

kcr fA 

We express the operator da in Hi as da = E^^i/ 
and /E = {^Ada\v)■ Here the superscript cr in de¬ 
notes that this SBO, when acting on a state, decreases 
the number of spin-cr electrons by 1. It is a Grassmann 
odd operator. We use without the superscript for 
general SBOs with unspecified quantum numbers. For 
Grassmann even (odd) Aap, its commutator (anticom¬ 
mutator) with Cka is zero. The SBO formalism has been 
used in the study of the Heisenberg model with large 
spini^ Recently this formalism is employed to combine 
the GF EOM truncation approximation with the exact 
diagonalization method to develop a new impurity solver 
for the dynamical mean-field theoryj^ 


A. Bare Strong-Coupling Expansion to Vj? Order 

Below, we first carry out the bare EOM series expan¬ 
sion for G(A'^p\A^\),^. The local single particle GF is 


fll' 


{A'^p, da>} — ^2 


(52) 


The coefficients and N^p^^^ read 

= dfiafvjS* + dvpfai^ i 
Hali,iLU — diiafpn + ^a' 

Using these definitions, we have 

[Alp,Hi]=Y,Y.^,a 


(53) 






r' A 


ka' fAV 


(54) 


obtained as 

GMap\A% 

_ ^i “1“ ^7a{'^(5/3)z 

ijJ + Ea — E^ 

kcr' fAiA H 

T/ 


ka' fiiA 


For i = 1, the above equation involves new GFs 
Go (^A^,,Cka'\A’22j and Go • They are 

zero because the impurity and bath are decoupled in Hq- 
Using the sum rule = (1)* = <^ 4 . 0 . the self- 

consistent solution gives {Aa 0 }i = 0. Therefore we get 

(-’l{A'^ap\-d2\)‘^ — 0 - 


(56) 
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In general, Gi{A'^p\A'^\)^ = 0 for i odd. This is because, 
starting from the impurity, an electron has to hop an 
even number of times to come back to the impurity and 
a local propagator contains only even powers of 14 • 

For J = 2, Eg. dS^ has the unknown averages of the 
type {Aai 3)2 and new first-order GFs Gi{A^i,Ckc,'\A‘'^\)^ 
and Gi{c\^,A^^\A'^^)^. The former is to be solved self- 

consistently with G 2 {Af^y\A^\)i^ later. For the latter, we 
can write down their EOMs and solve them similarly. 
Note that here is a Grassmann-even operator. In this 
process, new zeroth-order GEs will be generated. They 
can be expressed by the known quantities 

and (cj.g.Cfccr)o; using the fact that the impurity and bath 
are decoupled in Hq. The final results are 


^ Eu Cka' 


^ Vka' 
+Go{A-ys\A^^)t^ - 


^ T E^ Ejy (^ktr' 


(57) 


and 


W -|- E^ — Ey -|- Ckcr' 


— G[){A^s\A\^^)i^ 


Vk.' 


E L'^y^.ySi^k<7')0 

w -f E^ — Ey + Cka' 


(58) 


In the above equations, the newly introduced coefficients 
L^y^xr and are defined as 

\Afiy, dl^] = ^ ( L^y ^^^Axr', 

At 

[Afj,y,da] = 'y^Hyy^^AxT- (59) 

At 

Their expressions are 

Gfiy^Xr = ^kt^frv ~ ^urffiX', 

H;^,Xr = Sf^xf^r - (60) 

The four averages {A^^Ckc,)i, {AsyCko)i, ( 4 <t^m 7 )i 
and {(^f.^Asv)i in Eg. (1571) and Eg. (1551) need to be cal¬ 
culated from the GFs like Gi{cka\Afj,^)yj etc. using the 
fluctuation-dissipation theorem. For (A^^Cfeo-)!, we solve 
the EOM for Gi{cka\A^^)yj and get 


Gl (Cka I dl^7 ) u) 

^k<7 f -f + O/i) 

E^ — Efi -|- Cfecr 


I 

W - efecr 


1 

UJ E-y - Efi 

( 61 ) 


which gives 


{Af^^Cka') 1 

+ O/i) 

E^ Ef^ -\- Cka 


1 

+ 1 


1 

(62) 


Using {c].^Ayi~y)i = {A-yyj_Cka)i and the replacement (/i ^ 
i5 ,7 —>■ i/), the other three averages can be obtained. To 
carry out the A:-summations in Ea. (l55L it useful to in¬ 
troduce the following intermediate quantity, 


Kpi.^) = E 

k 


^ka {Aaf^Cka ) 1 
W - Cka 


/^a(ac + a/3) 
UJ -|- Ej^ — Ea 




(63) 


with 


,_<T _ A C A ,'17 Z7 A ^<t{uj) —r„{Ea — Ep) 

EPai^) — ^a{uj) A„{Ea Ep) ^PiEo^-Ep) ^ 

(64) 

Here, the two involved functions are 


r.(w) 


Aa{uj) 



A(e) 

UJ — c 


dc] 


A(e) 1 

w — e -|- 1 


de. 


(65) 


Likewise, the Hermitian conjugate of Eo. (l63l) is 


E 


^ka{(^k„Aap)i 

UJ - Cka 


®?a(w). 


( 66 ) 


With these preparations, we put Eg. (1571) and (1551) into 
Eg. (1551) . simplify all the terms and obtain 


G2(^a,al^7l)i7 

Sps^A 

a7 )2 H“ *^7a(-^(5/3)2 
UJ + Ea — Ej3 

(^uj + Eq, — -£//?) (clj + Ej — E§) 


(67) 


with given by 


^0.0,-^5 (^) 

= — E ^aP.usf^ui^fJ. + (ks)<P^y{uJ + Eyj, — Es) 

+ E ^^aP,-ivfvs{'^v + ‘^s)‘PysiuJ + Ej — Ey) 

vg' 

~ E ^ap,ki,sf^'yi^fi + — Efi+ Es) 

flG' 

+ E ^aP,yyfsu{‘^i' + ^s)^Sy{~‘^ — E^ + Ey). 

vg' 

( 68 ) 
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= (a^ + as) ^ ^a 0 ,'rufv 5 ^<y' {^ + — E^) 

vg' 

— (a^ + as) ^ -^a^,/i 5 //l 7 *ro-'(—w — E^ + Es) 

flG' 

— {a^ + as) ^2 — E,^) 

g' 

— (aj + as) ^2 ^2 — E^ + Ey). 

g' 

(69) 


Eas. (|67l) - (|69]) are the results for where 

Jafs-ysi^) involves and comes from the averages 

{(2.^Aap)i. Note that in the calculation we have made 
use of the Grassmann-odd properties of A'^^. Therefore 

these equations do not apply to GF G{Aai 3 \A 2 )^ with 
arbitrary Aap. 

Below, we focus on the particle-hole symmetric case 
and simplify these equations for numerical calculation. 
For the single impurity Anderson model considered in 
this paper, the four eigen states of the local impurity 
Hamiltonian hg = ^rio- are |1) = d||0), |2) = 

dJlO), |3) = |0) , and |4) = d|dj|0). The corresponding 
eigen energies are Ei = E 2 = €d, E^ = 0 , and E 4 = 
U + 2ed- There are total 16 SBOs Aap (a,/3 = 1 ~ 4). 
The impurity electron annihilation operator is expanded 
as d'f = A 31 + A 24 and dj, = A 32 — A 14 , meaning f 2 = 
f2 = 1 and /32 = —/i 4 = 1 and others are zero. In terms 
of electron operators. An = n-|-(l —nj^), A 22 = ( 1 —n-i-)n 4 ,, 
^33 = (1 “ n,-|-)(l — ni), and A 44 = n 4 ,n-f-. Their averages 
play important role in the T-dependence of GFs. 

The particle-hole symmetry implies the following con¬ 
ditions 


ed = -U/ 2 ; 

A^i-to) = A^iuj). (70) 


We have Ei = E 2 = —U /2 and Fg = £"4 = 0 to simplify 
Eos. ((501) and ((071) - (1001) . The local GF up to order 
reads G{da\dl)^ Ki Go{d^\dl)uj + G 2 {d^\dl)^, with 


GMdDu, 

G2{dM)u. 


1/2 ^ 1/2 

uj + U/2 ^ W-C//2’ 

Wr-1/2 W^-II2 W^jw) 
uj + U/2 lo-U/2 (w -b 17/2)2 
H7f(a;) W^iu) 

(w - [7/2)2 + (^ + U/2){uj - U/2) ■ 

( 71 ) 


The weights W))' (i = 1 ~ 5) for spin up are 


wl = 


IFJ — — (^22)2 + (^44)2; 


- -b (Aii )2 -b (^ 33 ) 2 ; 


= 2 -‘d52(a^) 




= 2 [ 7 ^ 14 “7^32(-a^) 


+ 2 [rt(‘a) + 04.(w) - r4,(-w) - A4 ,(w) + A4.(-w)]; 


1 r 


W^5'(w) = O E+'fi 2 {^) 


2 1 

+ ^ [-r.^(w) +r4.(-a;)]. 


(72) 


For spin down. 


wf = 
if/ = 


- + (^22)2 + (^33)2; 
2 + (^11)2 + (^44)2; 


W^iuj) 


2 1 


<pL(-^) - v>li{^) 


+ 2 + A-|-(w) - A^(-a;)]; 


= 2 [‘^l4(a^)-<d3i(-‘*^) 


+ 2 ^t(^) “ r-|-(—w) — A-|.(w) -b A-|-(—w)]; 


Wliu;) = - 


[rt(w) -rt(-w)]. 


(73) 


Note that only the averages of diagonal SBOs (Aaa) are 
involved in G(dcrjdDaj- This is due to the fact that both 
d^ = A 31 + A 24 and dj, = A 32 — A 14 consist of SBOs with 
nonoverlap subscripts. It is found that the relation H/f + 
IFf + W§ = rcr(w) holds here. The second condition 
in Eg. dTOl) implies rcr(—w) = —rg.(a;) and A„{—uj) = 
A 5 (a;) — rg.(tLi) which can further simplify the expressions 
for W:^{uj) (i = 3 ^ 5). Especially, we find IF^(a;) = 
IF 4 ^(w). 

Now let us consider the self-consistent determination 
of the averages {Aaa )2 (a = 1 ~ 4) which appear in IFf 
and IFf. They can be calculated from 02 (^ 311 ^ 3 ^)^, 

^ 2 (^ 321 ^ 32 ) 0 ;, 02 (^ 241 ^ 24 ) 0 ;, and G 2 (Ai 4 |A| 4 )o,, re¬ 
spectively. Using IFf (i = 1 ^ 5) dehned above, these 
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GFs read 


G'2(^3i|^3i)w = 

G2(^24|^24)‘^ = 

02(^321^32)'^ = 

G2(Ai4\A\^)^ = 


Wl-112 


7+7(4 

uj + U/2 


(w + [7/2)2 

7+2^ - 1/2 


7+44 

cq 

A 

1 

3 


(w - [7/2)2 

- 1/2 


7+7(0.) 

UJ + UI2 


( 0 ; + [7/2)2 

- 1/2 


7 + 40 .) 

UJ - [7/2 


(a; - [7/2)2 


(74) 


From these GFs, the self-consistent equations for {Aaa }2 
(a = 1 ~ 4) are obtained as 


(Aii)2 

(7133)2 

ePUi2 + 1 

g-/ 3 C //2 + 1 

(^44)2 

(^22)2 

g-/ 3 C //2 _|_ 

g/ 3 C //2 _|_ 1 

(^22)2 

(^33)2 

gdU/2 _j_ 1 

g-/ 3 C //2 _|_ 1 

(^44)2 

(Aii)2 

g-/ 3 C //2 _|_ 1 

g 0 U /2 _|_ 1 


/ wUlo) \ 
\{u2 + U/2YI 

/ Wjjco) \ 
\ (a.-[7/2)2/ 

/ \ 
\ (a.+ [7/2)2/ 

/ Wjicv) \ 
\ (a.-[7/2)2/ 


■ (75) 


In the above equations, the symbol {g{uj)) is defined as 
{g{uj)) = —1 /tt J^^lmg{uj + ig)l/{e^‘^ + l)duj. Note that 
the four equations in Ea. (l75ll are not independent. One 
needs to supplement the forth independent equation, i.e.. 


(7lii)2 + (^ 22)2 + (^ 33)2 + (^ 44)2 — ( 1)2 — 0. (76) 


These equations have to be solved numerically except for 
the case of paramagnetic bath where analytical solution 
is given below. In the limit T = 0, care must be taken in 
evaluating the right-hand side of Eq. ( 1771 ) and solving the 
equations, because the diverging /3 factors in the averages 
{Aaa )2 (a = 1 4) can make the numerical process 

unstable. In our numerical calculation, we isolate the 
most singular item from the above equations to stabilize 
the solution. 

In the paramagnetic phase, A-|-(a;) = Aj,(w). One finds 
lp\A—uj) = ipiy(uj) and cjL)- w) = Eos. ( 1771 ) and 

dZl give IFf = = 1/2, W^{uj) = Wi{uj) = F(w), 

and lF^(a;) = —F(a;). The self-consistent calculation 
of Aaa in Eg. ( 175 )) is no longer necessary. We obtain 
G(d,|4)^ = Go(d,|4)<.. + G2(d.|4)<.. + - and 


Go(d.|4)<.. 

G2(d.|4)<.. 


1/2 ^ 1/2 

uj + U/2 w-[7/2’ 

r(a;) F(a.) 

(w + G/2)2 (uj- [7/2)2 

(uj + U/2)(uj-[//2)' 

(77) 


This recovers the GF obtained by direct expansion in 
Ref. m. 


B. SE-resummation and CF-resummation 

In this subsection, we do the resummation for the 
second-order strong-coupling expansion obtained above. 
The most frequently used resummation method is via SE 
Eg. (175)1 . Although we have argued that SE resummation 
cannot solve the causality or the zero temperature diver¬ 
gence problem in general, here we will show the data for 
comparison. The second method that we will use is the 
CF resummation method Eq. (l20l) which is guaranteed to 
be causal but may have the zero temperature divergence 
problem. 

For SE resummation, from Eq. dm we obtain up to 
order 

= [7/2 + ([7/2)Vw; 

Sicr(w) = 0; 

E24cu) = Go\dM)Md,\di)u. - F4w). (78) 

The SE-resummed GF is then obtained by insert- 
ing this SE into the Dyson equation G {da\dl.)ui = 
[0//^(w) - E<^(w)] \ One gets2^ 

GsEida\di)ui 

^ uj + U/2- T^iiv) - Eo.(a;) - E2<.(a;) ' 

To do the CF resummation, we first expand Eq. dm 
into Taylor series of l/w. Note that the w-dependence in 
W[{uj) {i = 3,4,5) arises solely from the hybridization 
function A(e) which is proportional to V^. We therefore 
only expand uj in the denominators of Eq. dm and treat 
1+7 (w) (* = 3,4, 5) as constants. The obtained series is 
compared with the same expansion of GcF{da\dX)uj in 
Eg. (1701) . By requiring that for every n G [I,cx)], 0 ;“" 
terms in the two GFs agree on the level of V^, we obtain 
the expansion of the coefficients ao,ai,... and bi,b 2 ,---- 
It turns out that only two levels of fraction are sufficient 
to match the GFs up to V^. We get 


GcF{da\dX)u} 


Oo 


UJ + bi 


5 

ai 

UJ + b2 


(80) 


with coefficients 
Oo = 1; 

ai = ([7/2)2 ; 

b 4 = ^{W^ -W^)-T,{uj)- 

b 2 = - 7+2") - r.(w) + 21+5"(w). ( 81 ) 

Here, oq > 0 and oi > 0 guarantees that GcF{dcj\d\,)i^ 
has real simple poles and is causal. If we also expand uj 
in W[ {i = 3,4,5) to make the comparison, due to the 
continuous bath degrees of freedom, a CF with infinite 
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number of levels will be required to match Eq. (TflT) to 
order. 

Wi and W 2 enters the pole position of GsE(dcrldl)uj 
and GcridirldDu;, either through E 2 ( 7 (w) in Eg. (1791) or 
through 61 and &2 in Ea. dSOl) . Since they contain a /3 
factor, both GFs have the zero-temperature divergence 
problem as shown below. 

For the paramagnetic bath, SE-resummation and CF- 
resummation give respectively 


GsE(daldl)uj 


1 


3r(tj) 


(82) 


and 

GcF(daldl.)ui 


1 

a;-r(a;)-(|)V[c^-3r(cc)]' 


(83) 


Eg. (1821) was obtained by a direct expansion method in 
Ref. l28| and Eq. (IMl) was proposed there as an ad hoc rem¬ 
edy of the causality problem. It is seen that GsE{da\d}^)uj 
has a complex pole, leading to violation of sum rule 
in the local spectral function as shown below. 

GcF{dcr\dl.)ui has real simple poles only and it conserves 
the rum rule. Because of VEf = W 2 = 1/2 in the para¬ 
magnetic bath, in both results, the zero temperature di¬ 
vergence problem does not appear. 


C. Self-Consistent Strong-Coupling Expansion to 
Efc Order 


In this subsection, to remove the zero-temperature di¬ 
vergence problem in the bare expansion both with SE 
and CF resummation, we calculate G[da\d\)uj to order 
V)? using the self-consistent EOM expansion method. 

We split the Anderson impurity model HAim — Hq + 
Hi as in the bare strong-coupling expansion and use the 
same SBO definition. In the following, we use the self- 
consistent EOM expansion method described by Eq. m 
and (I24|) . The commutators and anti-commutators in¬ 
volved in the calculation are same as in Sec.IV.A. We 
will skip the calculation details whenever they are the 
same as before. 

The zeroth-order GF can be obtained easily from its 
EOM and we obtain 


Go{a-^p\a:^Ii 


{Asp) + 5ps{A 07 ) 

uj + Ed — Ej^ 


(84) 


Here {Aap) is with respect to the full Hamiltonian HAim- 
The first-order renormalized contribution is similar to 
Eg. (1551) with i = \ but without the average terms. 


Gi{A1^p\A''^\) 

= EE 


-Co 


kcr' 


OJ Ea — Ep 


-EE 

ka' fj-y 


Vka'Ki,, r 
UI + Ea — Ep 


(^Afj_,yCka'\A'^l'^ ^ 

( 85 ) 


The new zeroth-order GFs Gq and 

^0 can be solved from their EOM as 

ri ( A I /lO’t/ ~dvs{A'A^Ckai) + dj^{A \cka') 

Go [Af,,Cka'\A^A = -—-- -, 

\ / OJ U) Ck(7' “r En Ey 


( 86 ) 


and 


^0 ^ 


w -|- Cka' + En — Ey 


(87) 

Differing from the bare EOM expansion method, here the 
averages like {A'A^Cka') are with respect to the full Hamil¬ 
tonian HAim and non-zero in general. Putting Eas. (l 86 l) - 
(1571) into Eg. (1551) . we obtain the renormalized first order 
GF as 


Gi{Aap\A'^\)^ 


= EE- 

ka' fiy 


{uj + Ea 


-EE 

ka' k^y 


^ka' 

[uj + Ed 


Svs{A^j^Cka') + djfj.{A'^lcka') 
EjS^i^AJ ^ka' “t“ E^ Ey^ 

dys{(^ka'-^'^'i) — 

- Ep){u! -b Cka' + Efj_ — Ey) 

( 88 ) 


It contains contributions in all orders of Vka through the 
averages like {A^j^Cka')- 

The self-consistent EOM for the second order GF is 
solved in a similar way as the bare one and we obtain 


Mrp[i..Go ([A^.Cka ,, Hi] I^ 




= 2.2. 

ka' fJ-y 


(cJ “h Ed Ej^'jl^UJ €.ka' E E^ Ey') 


-EE 

ka' fiy 


Vka'KE^Go ([cL Hi]\A; 1)^ 

{ui + Ea — Ep){LO -b Cka' + E^ — E^) 

(89) 


Two new GFs of zeroth order appear and their EOM can 
be solved to give 


Go(K.Cfc,qi7i]|A;j)^ 
^ka' / ^ Jyr 

T 

+T.Y. 


pa" At 

-EE 

pa" At 


U) + Efx — Er 

^ p.y,XT {(*^T(5^A7 ^A7-^5r) Cpfj" ) 

^ H“ Ex Ej- ^pa" ^ka' 

-b SxjAst) (^^^iiCka'^ 


(a) Ex Ej- Cpcr" ^ka' 


( 90 ) 
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and 


proximation which is valid at the order 


_ _Y^ ^ jo'* ^v6{At^) + S^rjAsv) 


-EE 

per" At 

+EE 

per" At 


uj + E't — Ej^ 

^ p.u,\r {^rS-^X-f + (5A7^5r)^ 

+ Ex — Er + Cka' — ^per" 

^P<^" ('^Icr'E" + <5A7^5r)^ 


Lu + Ex — Ej- + efccr' + £’ 


per” 


(91) 


To complete the full self-consistent EOM series expan¬ 
sion, the averages on the right-hand sides of Eas. dS^ . 
dMl), dSOl) and dHI) are to be calculated from the corre¬ 
sponding full GFs, which themselves should be obtained 
from a finite order EOM expansion and CF resumma¬ 
tion. In this way, each average is obtained from a causal 
GF and contains /3 only on the exponent. This process 
has certain degrees of variance because one can choose 
the expansion order for the GFs used to calculate these 
averages. Given that the target GF G(A'^i^\A'^l)^ is pro¬ 
duced rigorously up to V^, different ways of calculating 
the averages amount to different approximations for the 
higher order contributions. 

Here, to avoid further complication, we will calculate 
the averages in the simplest way that keeps \A^\)ijj 

exact up to V^. For this purpose, the averages {Asp) 
and {^ 0 ,..^) in Eg. ((Ml) should be calculated from the CF- 
resummed GF GcF{A'^p\A’^\)^. Since the averages in 

^XAap\A’^^\)ui of Eg. (155)) already has a factor Vk in front 
of them, they will be calculated accurately to the 14 
level. For example, {A'A^Cka') can be calculated from 
the approximate GF G{cka'\A^'l)u: ~ Go{cko'\A'A^)^ + 
Gi{cka'\A'^jj))uj. Using the self-consistent EOM expan¬ 
sion, we obtain 

Goickcr'lA'Ajui = 0 ; 

GXcka'lA^l).. = -E^J2Ef3GoiA-o!p\A^^l)^. 

(92) 

In the second equation above, Go{A'^^\A'^j^)i^ is the 
renormalized zeroth-order contribution Eg. (1841) . Note 
that these contributions have real simple poles al¬ 
ready and the CF-resummation is not necessary. After 
{A'A^Cka') is obtained from the above GF, other averages 
in Eo. (|M)) can be obtained by subscript exchange, such 
as = {A^^Cka')*- 

The averages in G 2 (A^^|A^J)^ of Eqs.dlOl) and (IMI) 
appear on the level of V^. We will use a truncation ap- 


{AxjCpfjff Cka') SS {Ax^) {Cp(jH Cka’) Q b; 

(Cpcr"Cfccr'A a^) ~ {(^pa"^k(j')o{Ax"i) ■ (93) 

Similar decoupling approximations are used for other av¬ 
erages in Eo. dTO) and those in Eq. m- Putting these 
approximations into Eqs.((90]), (1911) . and Eo. (l89l) . the 
second-order GF is obtained as 




ra' 


UJ + Ea — E^ 


■ - u^ + E^-Ep - GMxAA^sV 


Xrer' 


- EE - u^ + E^-Ep - Go{A^AA^A^ 


pv Ter' 


Y^ Y^ ^aP,kyG%,XT^‘T'{ W Ep + EA .|. 

■ EE- co + E^-Ep -Go(Aa.|A^,).. 


pi' Xrer' 


(94) 


At this stage, we will make a further approximation 
to Gi{Aap\A'Jl‘Aui [Eq.([n21)] and G2(Aa;3 [Eq.([Ml)]. 
In these equations, Go(Aq,^|A^J);^ appears on the level 
of and we can make simplifications which are exact 
at 14 = 0, 


Go{Aap\Al^\)i^ 


^0(7 {Asp) + 5ps{Aaj) 

UJ -|- Ea — Ep 

{Aaa) + {App) 


Sa'ySpS 


■ Er, — Ea 


(95) 


That is, among the contributions higher than 14^, we ne¬ 
glect the processes that couple different SBOs and only 
consider the GFs that are diagonal on the BSO basis. 
With this approximation, Egs. dM)) and (IMl) are simpli¬ 
fied greatly. Putting Eo. d^ into Eo. dMl) and introduce 
the intermediate quantity similarly as in Eq. (1631) , one ob¬ 
tains 


= 


^ ^ (-^a/3Cfc, 

k 


- eka 
fpa {{Aaa) + {App)) 


UJ ■ 


Ea — En 


ifpaiuj). (96) 


is given by Eo. dM)) . Comparing the obtained 
Gi{A'^p\A^^Jju: with the bare EOM expansion result 
Eqs. d^-dMl), we obtain 


Gi(a-^|a;J), 


_ JaP,js{^) _ 

(cc -|- Ea — Ep) (uj -|- Ej — Es) 

( 97 ) 
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To simplify we put Eq.dM]) into 

Ea. dMl) and obtain 




_ ^aP,-is{^) _ 

(w + Ea — Ep) (oj + Ej — E^) 

(98) 


In the above two equations, ^a /3 ■ys(‘^) 

share the expression of JZp,-ys{^) (Eq-dMl)) and E^^^^g(uj) 
CEg. dMl) !. but with the substitution Ua —>• (Aaa) (a = 
1^4). Note that the renormalized zeroth-order GF 
Go{Actp\A^\)i^ in Ea. (l84l) is not changed and the aver¬ 
ages there need to be calculated self-consistently. 


Their CF-resummations read 


Gcf{A\^\AW)i^ — 

GcF{j^2i\^2i)t^ = 

Gcf(42|4I)- = 

GcF{A\^\Af^)u. = 


The quantities W[ {i = 1 
are given by 


{u + U/2)Wl -Wliuj)' 
(uj-U/2)W^ -Wj(uj)’ 

{u} + U/2)W^ -W^{uj)' 

-,, noi) 

(w - U/2) - ITi(w) 

~ 5) in the above equations 


Summing up Eqs. dH, (EZl), and (IMl) . we obtain one 
of the simplest self-consistent schemes of EOM expansion 
for G{A’^p\A'^l)^ which is exact to V^, 


G{AZ^\A% = 


{Asp) + SpsjA O 
(jJ + Efy - Er 


+ 


JaP,'ys{^) + ^aP^-ysi^) 


cx^^-yS V 


{ui Ea — Ep) (w + E^ — Es) 

(99) 


The averages of the type (A^p) in the above equation 
needs to be solved self-consistently from the full GF. Due 
to the same reason as discussed in the bare expansion, 
only averages of diagonal SBOs {A^a) (a = 1 ^ 4) are 
involved in G{dcr\d\)t^ . The advantage of the present 
self-consistent scheme is that it keeps the form of the bare 
GF expansion but only renormalizes the average values 
of the diagonal SBOs Aaa (a = 1 ~ 4). 

Let us now consider the self-consistent calculation of 
the remaining averages {Aaa) {a = 1 ^ A). These 
averages need to be calculated from the GF-resummed 
GFs GcF(^3i|^3i)tj, Gcf(^24|^24)<^: Gcf(^ 32|^32)<^ 
and Gcf(^i 4 |^i 4 )w Under the particle-hole symmetry 
condition Eg. (1701) . Eg. dMl) gives these GFs before resum¬ 
mation as 


g(4i| 4I).= 
G(44l4i)^ = 
G(42l4J)- = 

G(44i4i)^ = 


wl 


4h 

uj + U/2 


{uj + 17/2)2 

wj 


Wjiio) 

UJ-UI2 


{uj - 17/2)2 



4(02) 

uj + U/2 


{uj + 17/2)2 

4^ 


Wj;{uj) 


UJ-UI2 (w-17/2)2 ■ 


( 100 ) 


~ t 
lUi 

~ "t 

IU 2 

#4'^(a;) 

lU5^(w) 


and 

it/ 

1 ^ 2 '^ 

ITs/w) 

W4^{UJ) 

ITs/w) 


hs] 

hi] 

hi'^u{-^) - h3vi2{^) 

+7i 3 [rt(w) -I- A_[,(w) — A4_(—w)]; 

hi^u{u^) - h3^i2{-^) 

Ehi [r't-(w) -|- r|(a;) — r4,(—w) — A_|,(a;) -I- Aj,(—w)]; 


-hi 

(p\h^) + (p\^{-Uj) 

Ehs 

Pi2{^) + P32{-^) 

+hi 

-rb(w) -br|(-a;) -b 

+7i3 

1 

> 

■(- 

+ 

> 

■(- 

1 


( 102 ) 


hs] 

hi] 

hi^P^i-^) - hsplh^^) 

+h3 [r4.(w) + A-|-(a;) — A-|-(—w)]; 

hiplh^^) - hsplh-i^) 

+7i 4 [r4,(w) -b r-|'(a;) — r-|-(—w) — A-|-(a;) -I- A-|-(—w)]; 

hi plhuj) + 

-In 4i(w)+4i(-w) 

-I-7i 4 [r'|-(w) — r'|'(— uj) — A'|-(w) -|- A-|'(—w)] 

+^23 [A't-(w) — A'|'(—w)]. 

(103) 


Here, hp = {Aaa) + {App). The relation tT^(a;) + 
TTf(a;) + Wg{uj) = rcr(w) still holds. Gompared to 
the same quantities in the bare expansion Eg. (1721) 

and Eg. (1751) , the 1/2 factor there is replaced with the 
averages lap in the renormalized expansion. The self- 
consistent equations for the averages {Aaa) (a = 1 ~ 4) 
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are completed by 


^ r°o _ 2^ 

{■^Pp) — / lTOG'cF(^a/3l^Q^)‘^+i») 

^ J —oo ^ “r i 

(104) 

and J2p{^pp) = 1 - 

With the averages (Aaa) (a = 1 ^ 4) obtained, the 
single particle GFs G{d^\d\)u = G{AI^+ 

and = G{Ai^-Ai^\A^^-A[\)^ are calculated 

as 


G(d,,|4)^ 


Wf 

u + U/2 
, W^uj) 
(w - U/2) 


uj-U/2 {oj + 17/2)2 

W^iuj) 

2 {uj + U/2){uj-U/2)' 
(105) 


The CF resummation is then carried out to it in the same 
way as for the bare GF expansion Ea. dMl) and The 
result is denoted as Gsc{da\d^a)i^ '^6 obtain 


Gsc{dM)^ 


0-0 


cj + bi 


Oi 

a; + 62 


(106) 


with coefficients 

Oo = 1; 

01 = (U/ 2 f ; 

bi = j(w,^-W2n-rAco); 

b2 = - m) - ra(o^) + (107) 

To obtain these results, we have made use of the fact 
that — IT 4 ' oc and neglected them when compar¬ 
ing the l/w expansion of G{da\dl.)ui and Gsc{d<T\dl)uj- 
For the paramagnetic bath, Eos. (11061) and (11071) reduce 
to Eas. d^ and (ISTI) of the bare EOM expansion with 
CF-resummation. Note that the inverse order: First do¬ 
ing CF-resummation for G{Aap\A'^l)t^ and then sum¬ 
ming them up to produce G(dcr|dJ-)<^i does not work. 
This is because some components, e.g., G{A 3 i\Al^)t^ and 
G(A 24 |^ 3 i)cj, starts from l/w^ in the l/to expansion and 
the form of CF Eg. (11061) does not apply. 


D. Numerical Results 

In this subsection, we present numerical results for the 
formula obtained in previous subsections. We compare 
results obtained from the three different combinations of 
second-order strong-coupling expansions and resumma¬ 
tion methods: bare EOM expansion with SE resumma¬ 
tion (bSE), bare EOM expansion with CF resummation 
(bCF), and self-consistent EOM expansion with CF re¬ 
summation (SC). All these results are compared with the 



FIG. 1: (Color online) The impurity density of states 
for various U’s, obtained using (a) bSE; (b) bCF; (c) SC; 
and (d) NRG. From top to bottom at small the oj regime: 
U = 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, respectively. Other model pa¬ 
rameters are F = 0.1, Acj = 0.0, ed = —U/2, and T = 0.1. 
NRG parameters are A = 3.0, Ms = [256, 280], log-Gaussian 
broadening parameter B = 0.08, and Nz = 8 . 


numerical renormalization group (NRG) data, which is 
believed to be accurate at the low and small-frequency 
regimes. We use a Lorentzian hybridization function for 
the Anderson impurity model. 


Acr(a;) 


Fw 


2 

c 


{oj -I- (TAa;)2 -I- w/ ’ 


(108) 


Here, Wc = 1.0 is the energy unit. F is the hybridiza¬ 
tion strength. The spin-dependent energy shift Aw is in¬ 
troduced to simulate the spin-dependent bath energies, 
cr = -1-1 for spin up and cr = — 1 for spin down. Aw = 0 
gives paramagnetic bath while Aw yf 0 mimics the spin- 
polarized bath. In this paper, we only study the particle- 
hole symmetric point Ed = —U/2. 

NRG calculation of the local density of states (LDOS) 
for this Anderson impurity model is done with the full 
density matrix formalism,— supplemented with the SE 
trick— To discern the sharp Hubbard band in the large- 
U regime, we use a small broadening parameter and aver¬ 
age the LDOSs with = 8 interleaved discretizations— 
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FIG. 2: (Color online) (a) The sum rule of impurity density 
of states as functions of U in paramagnetic bath; (b) double 
occupancy as functions of U. The symbols with guiding lines 
are bSE (circles), bCF (up triangles), SC (squares), and NRG 
(down triangles). The model parameters are F = 0.1, Alj = 
0.0, ed = —1//2, and T — 0.1. NRG parameters are same as 
in Fig.l but without the z-average. 


The final LDOS from this standard procedure is believed 
to be accurate at least for low temperature and in small 
frequencies. The thermodynamical quantities {ria} and 
(n-^ni) are obtained from the respective spectral func¬ 
tion by frequency integration. 

Figure 1 presents LDOS at T = 0.1 and various C/’s for 
the paramagnetic bath Alu = 0. The data are for bSE, 
bCF, SC, and NRG in Figs. 1(a), 1(b), 1(c), and 1(d), re¬ 
spectively. Among the three combinations of expansion- 
resummation method, LDOS of bSE has an unphysical 
dip at cu = 0 for any non-zero U. bCF and SC give 
out identical LDOS because SC reduces to bCF in the 
particle-hole symmetric and paramagnetic situation. The 
height of the central peak obtained from bCF [Fig. 1(b)] 
and SC [Fig.1(c)] decreases gradually with increasing U. 
This behavior, being consistent with NRG, is correct for 
temperature higher than the Kondo temperature. In 
all the obtained data, the Hubbard peak positions are 
slightly larger than U/2 due to the hybridization shift. 
In the large U limit, all LDOSs tend to the atomic form 
Pati^) = l/2i5(a;—C//2)-|-l/2(5(a;-|-C//2) which is expected 
when U is much larger than bath band width. Compared 
to the NRG curve in Fig. 1(d), qualitative agreement is 
reached by bCF and SC in both the small- and large-C 
limits. The deviation is stronger in the small-w regime 
for intermediate U values 1.0 < U < 3.0. 

The good quality of LDOS at small U obtained from 


strong-coupling expansion is a consequence of resum¬ 
mation which effectively extends the validity range of 
the series expansion. Actually, all the three expansion- 
resummation schemes give exact GF at C = 0. It is 
observed that for large U values, the Hubbard peaks are 
significantly sharper than NRG results. The neglecting 
of higher-order contributions of hybridization may lead 
to sharper Hubbard peaks, but we believe that the main 
reason for this discrepancy is the poor energy resolution 
of NRG at high energies. It is known that NRG tends 
to over broaden high energy peaks. Indeed, the height 
of Hubbard peaks increases when we use smaller broad¬ 
ening parameter, larger and keep more states. Here 
we have used the log-Gaussian broadening. Sharp fea¬ 
tures have been obtained from Gaussian instead of log- 
Gaussian broadening and Nz up to 32 

Figure 2 presents the sum rule and double occupancy 
as functions of t/ at T = 0.1 and Aw = 0.0. In Fig.2(a), 
integration of the LDOS in Fig.l are compared. As ex¬ 
pected, LDOS from bSE does not obey the sum rule, 
while those from bCF and SC obey it at a precision 10“^. 
The tiny deviation is due to numerical error. NRG result 
fulfills the sum rule at machine precision. 

The double occupancy {n^n^) can be calculated in var¬ 
ious ways. For bSE and bCF, it can be calculated directly 
from the single particle GF as 

1 1 
'TT J —QQ ^ “r -L 

G{nada\d'l)uj = —G{da\dl)ui^a{u}). (109) 

One can use either a =t or a =), in the above equations. 
For SC, besides the above equation, one could also use 
= (A 44 ) since (A 44 ) has been obtained in the 
self-consistent calculation. Different ways of calculating 
have relative deviations smaller than 5%. In this 
paper, for bSE, bCF and NRG, we use Eg. (11091) with 
a =t and for SC we use (n-t-nj,) = (A 44 ). 

In Eig.2(b), the bSE result for the double occupancy is 
much smaller than the other three, and even slightly neg¬ 
ative near U =1.0. The results of bCE and SC agree with 
that of NRG at quantitative level. Similarly to LDOS, 
the agreement of the bCE and SC results with NRG is 
better in both the small- and large-t/ regimes. 

In Eig.3 and Eig.4, we focus on the Anderson impu¬ 
rity model with a spin-polarized bath Aw = 0 . 2 , all at 
an intermediate temperature T = 0.1 = F. In Fig.3, the 
[/-dependence of LDOS is shown. Since p\.(ui) = p^{—u}) 
is obeyed very well, here we only show P't-(w). All the 
LDOS curves have an asymmetric shape due to spin po¬ 
larization in the bath. Similar to Fig. 1(a), bSE curves 
have unphysical dips at w = 0 for nonzero U. the bCF 
and SC results are similar, but no longer identical in the 
case of magnetic bath. It is seen that the bCF result has 
more asymmetry in the upper and lower Hubbard bands 
than does the SC result. 

In Fig.4(a), the sum rule of LDOSs is analyzed. bSE 
has an incorrect sum rule while bCF and SC fulfill it 
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FIG. 3: (Color online) The impurity density of states 
for various t/’s in a spin polarized bath, obtained using (a) 
bSE; (b) bCF; (c) SC; and (d) NRG. From top to bottom at 
the small-a; regime; U = 0.0, 0.5,1.0, 2.0, 3.0, 4.0, respectively. 
Other model parameters are F = 0.1, Auj = 0.2, = —V12, 

and T = 0.1. NRG parameters are same as in Fig.l. 


perfectly. In Fig.4(b), the impurity electron occupancies 
(n-|-) (filled symbols) and (nj,) (empty symbols) are shown 
as functions of U. Due to the incorrect sum rule, bSE 
gives a total occupancy less than half-filling. In contrast, 
(n-|-) — 1/2 = — ((nj,) — 1/2) is preserved in the results 
of bCF, SC and NRG. Being consistent with the larger 
asymmetry in LDOS, bCF gives qualitatively larger mag¬ 
netization M = Kn-i-) —(nj,)! than SC. Compared to NRG 
data, the bCF result agrees better for U > 3.0 while the 
SC result has better behavior for U < 1.0 and all curves 
are non-monotonic. The SC result for M is appreciably 
smaller than NRG even in U = 5.0. 

This quantitative difference can be traced back to the 
approximation scheme Ea. dMl) and (IMl) used for the self- 
consistent calculation of averages. In those equations, 
the atomic-like truncation scheme weakens the influence 
of the asymmetric bath on the impurity and leads to 
smaller M. In other words, although the single particle 
GF being exact up to V^, the averages are actually eval¬ 
uated with respect to a ground state or a density matrix 
which is accurate to lower order. As will be detailed later, 
the significant error in the large-17 regime hints that not 
only T/U, but also F/T should be small to guarantee 



U 


FIG. 4: (Golor online) Physical qnantities as functions of U 
at fixed T = 0.1 for the Anderson impurity model with spin 
polarized bath, (a) Sum rule; (b) averages (n-|-) (filled sym¬ 
bols) and (nj,) (empty symbols); and (c) double occupancy. 
The symbols with guiding lines are bSE (circles), bGF (up 
triangles), SG (squares), and NRG (down triangles). Model 
parameters are F = 0.1, Auu — 0.2, td = —U12. NRG param¬ 
eters are same as Fig.l but without the a-average. 


quantitative accuracy in the present expansion scheme. 

In Fig.4(c), the double occupancies are shown as func¬ 
tions of U. They look similar to the paramagnetic case 
while the bCF result has some drawback in the magnetic 
case: (n-j-nj,) exceeds the upper limit 1/4 in the small U 
limit. This reflects that although the single particle GF 
obtained from bGF is exact at U = 0, the two-particle 
GF is not. When U is small, something is qualitatively 
wrong in the higher order GFs obtained from bGF. In 
contrast, SC result agrees with NRG much better from 
small to large U. 

To study the temperature dependence of these results, 
in Fig.5 we show LDOS on log-scale atU = 3.0, Aw = 0.2 
for various temperatures. In Fig.5(a), it is seen that 
LDOS from bSE are normal at T = 0.1 (black solid line) 
(the dip at w = 0 still present). As temperature low¬ 
ers further, the lower Hubbard band stays around —11/2 
while the upper Hubbard band begins to split and the 
higher branch moves to larger values. Figure 5(b) shows 
the LDOS from bGF. In the limit T = 0, both the lower 
and the upper Hubbard bands move towards ±oo, with 
enhanced weight transfer from the lower Hubbard band 
to the higher one. Numerically we find that the position 
of these peaks is proportional to 1 /T in the small T limit. 
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FIG. 5: (Color online) The impurity density of states p^{oj) 
for U = 3.0 and various T’s in a spin polarized bath, obtained 
using (a) bSE; (b) bCF; (c) SC; and (d) NRG. The curves are 
for T = 0.1 (solid line), T = 0.03 (dashed line), T = 0.01 
(dotted line), and T = 0.005 (dash-dot line). Other model 
parameters are F = 0.1, Ao; = 0.2, and ed = —U12. NRG 
parameters are same as in Fig.l. 


They are the unphysical features due to the (3 factors in 
the bare expansion of GF. After either bSE or bCF re- 
summation, these factors enter the denominator and in¬ 
fluence the position of poles. In the second-order bare 
EOM expansion Eq. dZID, it is Wi and FF^ that contain 
the /3 factor via {Aaa )2 {a = 1 ^4). Eor the param¬ 
agnetic bath, Wi = W 2 = 1/2, the /3-dependent terms 
cancel and the problem does not appear. For the mag¬ 
netic bath, we numerically find that (^ 11)2 and {^ 22)2 
are proportional to /3, while (^ 33)2 and (^ 44)2 are almost 
independent of /3. 

The LDOS from SC shown in Fig.5(c) has very weak 
temperature dependence. As T decreases from T — 0.1, 
the Hubbard peak positions move weakly and converge 
to ±1.8. The weight distribution does not change much 
down to T = 0.005. Compared to the NRG results in 
Fig.5(d), SC gives quantitatively correct peak position 
of Hubbard bands. What is missing in the SC results is 
the central Kondo resonance at small T and the continu¬ 
ous weight transfer from the lower Hubbard band to the 
higher one as T decreases. This is again a consequence 
of the atomic-like truncation scheme used in Eq. (1931) and 
Eq. dMI) . Overall, in the low temperature limit, bSE and 
bCE have diverging positions of Hubbard peaks, while SC 
maintains correct peak position but has weaker spectral 



FIG. 6: (Color online) Physical quantities as functions of T 
in log scale at 17 = 3.0 for the Anderson impurity model 
with spin polarized bath, (a) Sum rule; (b) averages (n^) 
(filled symbols) and (n^) (empty symbols); and (c) double 
occupancy. The symbols with guiding lines are bSE (circles), 
bCF (up triangles), SG (squares), and NRG (down triangles). 
Model parameters are F = 0.1, Acj = 0.2, ed = —U12. NRG 
parameters are same as in Fig.l but without the 2 -average. 


weight transfer compared to NRG. In the high tempera¬ 
ture regime T ^ T, the shape of LDOS from bCE and 
SC agree well with NRG. 

In Fig. 6 , the sum rule, electron occupation, and the 
double occupancy as functions of T are presented. The 
sum rule in Fig. 6 (a) shows that bSE result is incorrect for 
all temperatures, while bCF and SC results keep at unity, 
as expected. The electron occupations in Fig.5(b) show 
significant difference among the four results at low T. 
The result of bSE does not fulfill particle-hole symmetry 
and (nj,) exceeds unity. bCF and SC give qualitatively 
correct results. At the low temperature limit, bCF gives 
out a fully polarized impurity M = 1.0, consistent with 
the enhanced asymmetry in LDOS at the low-T regime. 
In contrast, SC gives much weaker polarization, with M 
saturating to 0.084 at T = 0.005, much smaller than 
the NRG value 0.742. At the high-temperature regime 
T bCF and SC results agree well with NRG. 

Double-occupancy results are shown in Fig. 6 (c). The 
result from bSE is much smaller than the others and 
slightly negative around T = 0.4. At low temperatures, 
bCF produces a qualitatively wrong result: {n-fUi) tends 
to negative for T < 0.02. In contrast, the SC result de¬ 
creases with lower T at high temperatures and reaches 
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FIG. 7: (Color online) (a) and (b): The differences between 
the impurity magnetization M = | — nj.)| obtained from 

expansions and NRG, as functions of T. The parameters 
are (a) U = 0.05 and (b) U = 3.0, at Atu = 0.2. The 
filled symbols are for bCF and the empty symbols are for 
SC. Squares, circles, and triangles correspond to F = 0.03, 
F = 0.1, and F = 0.3, respectively, (c) Checking Friedel sum 
rule at T = 0.005 and F = 0.1 for small U values. Filled 
symbols are A.t'(0)Pt(0) empty symbols are sin^(7rn.|-). 


a constant {n^n{) = 0.03 for T < 0.4, close to the NRG 
value 0.039 at the low-temperature limit. At high tem¬ 
peratures, the double occupancy from bCF and SC agree 
well with that of NRG. 


V. DISCUSSION AND SUMMARY 

We first discuss the validity range of the expansion 
schemes. The quality of an expansion is controlled by 
the small parameter A in H, which in our strong-coupling 
expansion for the Anderson impurity model is character¬ 
ized by the hybridization strength F. Roughly speaking, 
the validity range of the expansion is F much smaller 
than all the other energy scales, including the interac¬ 
tion U, temperature T, and frequency w. In Figs.l ~ 4, 
T = F = 0.1 is used which is on the boundary of the va¬ 
lidity range, and the results have significant errors. It is 
seen from Fig.6 that the agreement between the expan¬ 
sion and NRG is good for T F and C/ F. Figures 


3 and 4 shows that the agreement is also good for small 
U ^ 0 due to the C/^-order accuracy acquired from the 
CF-resummation. 

Figure 7(a) and 7(b) demonstrate that F indeed sets 
in as the breakdown temperature scale in the present 
strong-coupling expansion approaches. Taking the im¬ 
purity magnetization M = (n-|-) — (nj,) as an example, we 
show the difference between the expansion and NRG re¬ 
sults Mexp — Mmrg as functions of T /F, for U = 0.5 ~ F 
in Fig.7(a) and for U = 3.0 F in Fig.7(b). For F 
values ranging from 0.03 to 0.3, the magnetization from 
both bCF and SC approach the NRG values at T/F 1, 
and deviate significantly at T/F < 1. Similar behavior 
is observed in the double occupancy curve (not shown 
here). In Fig.7(c), we examine to what extent the Friedel 
sum rule /S.a{Q)pa{0) = sin^(7rncr) is satisfied at low 
temperaturesIt is seen that for both spin polarized and 
un-polarized cases, this relation is fulfilled exactly only at 
[7 = 0. Approximate fulfillment in seen in {7 < F = 0.1. 
Since the Friedel sum rule is a consequence of the Fermi 
liquid ground state, it is not expect to be fulfilled well by 
the strong-coupling expansion which starts from the lo¬ 
cal moment limit. However, due to the CF-resummation 
method used in our approach, GF acquires the correct U 
term and the Friedel sum rule is satisfied in the regime 
U <T. 

To summarize, for T ^ F, both bCF and SC schemes 
are quantitatively accurate in the regime U ^ T, and 
produce a smoothly interpolation between large U and 
small U. For T < F, accurate results can only be ex¬ 
pected for ^ 0. Since the Kondo scale Tr is much 
smaller than F for large U, the Kondo resonance cannot 
be described by the the present strong-coupling expan¬ 
sion method. 

Although the SC scheme improves over bSE and bCF 
on the causality and the zero-temperature divergence 
problems, the appearance of a breakdown scale F in 
the temperature axis is an obvious shortcoming, sim¬ 
ilar to the strong-coupling expansion for the Hubbard 
model i^^d^ Considering that for the Anderson impurity 
model, much better results are available from meth¬ 
ods such as NRG, QMC, and functional renormaliza¬ 
tion group method^, the present expansion schemes re¬ 
ceives only partial success at this stage. However, the 
results are amenable to further improvement. The ap¬ 
pearance of the breakdown scale in temperature and the 
zero-temperature divergence problem are related to the 
fact that 14 = 0 is a singular point in the ground state 
of H with a spin-polarized bath. One could avoid these 
problems by selecting a suitable Hq whose ground state 
is continuously connected to that of H. This issue will 
be studied in the future. 

Below, we discuss some distinct features of the present 
expansion method compared to existing theories. The 
present approach is universal in the sense that it has no 
requirement on the form of Hq. For any Hq that is ex¬ 
actly solvable, i.e., either its GF EOM closes naturally, 
or its eigen states and eigen values are obtainable, series 
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expansion of GF can be constructed in a unified frame¬ 
work. Hq- and i/i-specific diagrammatic rules are not 
needed in this method. By using the self-consistent EOM 
expansion supplemented with CF resummation, causality 
of GF is guaranteed and the zero-temperature divergence 
problem removed. The resulting GF has extended range 
of validity. Therefore, for those Hamiltonians that both 
Hq and Hi are exactly solvable, by expanding GF from 
the two limits and comparing the results, one can obtain 
reliable knowledge in both the small and large interac¬ 
tion regimes. In principle, expansion around a cluster 
or impurity Hamiltonian Hq is also possible. This could 
provide a possible alternative derivation of the cluster 
perturbation theory^ or dual fermion dynamical mean- 
field theory.— 

The present approach is distinctive in that arbitrary 
double time GF can be expanded in the same framework. 
In the traditional methods, calculation of GFs of more 
than one particle is a laborious task. Here, due to the 
universality of the formalism, multiple-particle GF can 
be expanded in a way parallel to one-particle GF. For 
an example, the strong-coupling expansion method used 
in this paper for calculating the single particle GF can 
also be used to produce the dynamical spin- or charge- 
correlation functions, with only slight modification in the 
formalism. 

Next, we discuss possible improvement and exten¬ 
sions of the present approach. As seen in Fig.4(b) 
and Fig.6(b), (ua) from the second-order self-consistent 
strong-coupling expansion deviates significantly from the 
NRG result, even for U as large as 5.0. Also, the 
temperature dependence of M is too weak. Obviously 
there is much room to improve the result. To calcu¬ 
late the averages of the type {AZ^ka') (in Ea. (l88l) l and 
{Ax^Cpai'Cka') [in Eas. dTOl) and (ITO ]. instead of using the 
atomic-like truncation scheme, we can carry out the self- 
consistent EOM expansion for GEs G{cka'\A^jj^)ui and 
G{cka' \A\^Cp„ii)^ to V)? order and calculate the averages 
from the CF-resummed GFs. Also, a suitable selection 
of Hq may help remedy the zero-temperature divergence 
problem and remove the breakdown scale in temperature. 

A by-product of the present method is the EOM for 
the n-th order residue Tn{A\B)t^, Ea. dTTl) or Ea. (E5|) . It 
could be employed to produce higher order modifications 
to the series up to Gn{A\B)^. The EOM of Tn{A\B)^ 
is formally similar to that of the full GE G{A\B)^, ex¬ 
cept that the lower order contributions have been sin¬ 
gled out. In principle, it can be solved approximately by 
standard truncation schemes. This provide possibilities 
of constructing new types of truncation approximations 
which are exact up to order A", or developing improved 
GF resummation formulas with a terminator— 

An ideal expansion scheme may be that the resulting 
GF is exact simultaneously to and V^U'^, and 

hence accurate in both the weak- and strong-coupling 
limits. As it is difficult to realize in traditional perturba¬ 
tion theories, it is apparently possible to achieve this goal 
in the EOM expansion method. One could first carry out 


the weak-coupling expansion to obtain G^‘^'> {da\d^fj)u], and 
then carry out the strong-coupling expansion to V)? order 
for the residue T 2 {da\d}^)uj by employing its EOM. The 
resulting GF G^^’^\d^\d^^)i^ satisfies the above require¬ 
ment. In practice, however, the resummation method 
suitable for such an expansion is yet to be developed. 

Direct multiple-variable expansion is also possible 
within the framework of EOM expansion. Eor an exam¬ 
ple, the splitting of Hamiltonian H — Hq + XHi + 9H2 
can be used to generate a GF expansion G{A\B)^ — 
yO^Gij{A\B)^. If we choose Hq = J 2 ka^kacl^Cka, 

Hi = J2ka'^ka{cl„da + d^Cka), and H 2 = Un^n^ + 
the expansion up to lowest several orders can 
be obtained with ease. A subsequent resummation can 
be used to produce a physically meaningful result, be¬ 
ing correct in both the strong- and weak-coupling limits. 
If the full Hamiltonian is treated as a perturbation, the 
self-consistent expansion will produce a moment expan¬ 
sion of GF, while the bare expansion is equivalent to the 
simultaneous moment and high temperature expansions. 

The same strategy of expanding the double time GF 
can be extended straightforwardly to other GFs, if only 
the EOM formalism also applies there. For an exam¬ 
ple, the Keldysh GF describing the non-equilibrium pro¬ 
cess can be described by EOM. The present EOM-based 
expansion method can be extended to calculate the the 
Keldysh GF. 

From our demonstrative calculation for the weak- as 
well as strong-coupling expansion, it is clear that the 
present method also has some shortcomings. For most 
of the models, it is difficult to obtain explicit expansion 
higher than second order because the complexity of cal¬ 
culation increases very fast with order. This feature is 
common in every expansion method such as the Feyn¬ 
man diagram for weak-coupling expansion, Metzner’s di¬ 
agram for strong-coupling expansion;^ or Dai’s direct 
expansion— In these techniques, the time ordering and 
multiple integrals will complicate the problem. With 
the aid of computer algebra, we hope that higher or¬ 
der GF could be obtained, similarly to the situation of 
strong-coupling expansion— Another drawback of the 
present approach is that the partition function can only 
be obtained indirectly by using the coupling constant in¬ 
tegral method. The calculation of free energy is impor¬ 
tant for studying thermodynamical properties and con¬ 
structing the conserving approximations. The present 
EOM-based expansion basically expands the excitation 
energies instead of the eigen energies. We have not yet 
found ways to construct the direct expansion of partition 
function. Finally, differing from the diagrammatic meth¬ 
ods where a diagram in arbitrary order can be evaluated 
directly, in the present method, series can be generated 
only recursively and calculated order by order. 

In summary, we have presented an EOM-based method 
for doing series expansion of double time GEs. We de¬ 
veloped both the bare expansion and the self-consistent 
expansion formula. Using this method, we carried out 
the second-order weak-coupling expansion as well as the 
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strong-coupling expansion of the single particle GF for 
the single impurity Anderson model. For the weak- 
coupling expansion, Yamada’s SE up to C/^ is obtained. 
For the strong-coupling expansion, we obtained results 
from three different expansion-resummation schemes: 
the bare expansion with SE resummation, bare expansion 
with CE resummation, and the self-consistent expansion 
with CE resummation. The latter overcomes both the 
causality problem and the zero-temperature divergence 
problem. We found that although they agree with NRG 
well in the large-G and T ^ T regimes, quantitative 
accuracy is not achieved at low temperature. Some fea¬ 
tures of this new approach and possible extensions are 
discussed. 
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The free LDOS is obtained as /Oocr(e) = 
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Appendix A: Calculation of Gi{n^dcr\dl^)u, in 
weak-coupling expansion 

In this appendix, we calculate Gi{ngdcr\dl.)ui using the 
bare EOM expansion. Straightforward calculation with 
the right-side EOM gives 


Gi(n^d,|4)^ = Go(d.|4)-x 

4“ UGQ(^Ti^dfj\'Eiad^^uj o^crGoIj . 

(Al) 

This equation involves a zeroth-order two-particle CE 
Go(n^do-14)i.j and a three-particle CE Go{nada\ng- 4).. 
Go(n<,d,,|dt),, can be solved easily by its right-hand side 
EOM as Go(n<fdcr|dJ.),^ = {na)oGo{da\dl)i^. Direct EOM 
for the three-particle CE Go{nada\ngd}^)uj will lead to 
new three-particle GEs and the closure of the hierarchy 
is slow. So we first diagonalize the unperturbed Hamilto¬ 
nian Hq lEg. d^Sl) in the main text) in the single-particle 
space and obtain 


Tdo — ^ Esa^\(j^sa ■ (-^^) 

S(7 

Here s is the single particle orbital index. We assume 
= Y.s^saascj with |dscrp = I- Using the quasi¬ 
particle GF of ido 

Go(a..|a4,)- = EEi, (A3) 

^ Ss<7 


where -^suv,s'u'v' — The EOlVI 

for Go{algau<jaya\al,^au'aal,^)ui gives 



^ T £S(T ^U<7 ^V(T 


(A 6 ) 


The nominator of this GF is easily calculated as 
Syyf Sus' ^su'‘fT'VO'(^‘^ud' ) T ^w'^us'^su'^sai^^ H“ 

^us^u's^'^sa'^s'a ■ Here Tlsa — (^s(7^S(7')o — T 

1). Putting it into Eq. (|A 6 p and (|A5p . one obtains 


_ \^sa\ \hua-\ l^val ^ \ \ ^ _ M 

— / \'^UO' ^sa) T ^S<7 (1 '^ug)\ 

SUV ^ ~ ~ 

+ (nCT)oGo(do'|d^)(^. (A-f) 


In terms of the free LDOS, this equation is written as 


GQ{nifda\n^dl)^ = (n^)oGo(d^|d^), 
P0s{El)P0s{E2)P0a{Ez) 


w -I- ei - £2 - 63 


F{ei,e2,e3)deide2de3, 

(AS) 


with F(ei, 62 , 63 ) = - riej -f (I - Fi¬ 

nally, putting Ea. dASIl into Ea. dAll) and using the func¬ 
tion Ka{uj) defined in Ea. (l38ll . we get 


Gi{ngda\d'l)uj 

= [{ns)i + UK„{uj)] Go(dcr|d^)a; 

\U{Tla}o Q^cr] Gq (d^r |d],.)^;. 

(A9) 


This completes the calculation of Gi(nffdcrldl)uj- 
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